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Abstract 

A class of implicit Total Variation Diminishing (TVD) type algorithms suitable for transonic and 
supersonic multidimensional Euler and Navier-Stokes equations has been extended to hypersonic 
computations. The improved conservative shock-capturing schemes are spatially second- and third- 
order, and are fully implicit. They can be first- or second-order accurate in time and are suitable 
for either steady or unsteady calculations. Enhancement of stability and convergence rate for hyper 
sonic flows is discussed. With the proper choice of the temporal discretization and suitable implicit 
linearization, these schemes are fairly efficient and accurate for very complex two-dimensional hy- 
personic inviscid and viscous shock interactions. This study is complimented by a variety of steady 
and unsteady viscous and inviscid hypersonic blunt-body flow computations. Due to the inherent 
stiffness of viscous flow problems, numerical experiments indicated that the convergence rate is in 
general slower for viscous flows than for inviscid steady flows. 


I. Motivation and Objective 

Most shock- capturing methods are either inefficient for practical computations or only valid for 
transonic or supersonic perfect gas calculations. For hypersonic, perfect gas, equilibrium real gases or 
nonequilibrium flows, improvement and modification to existing methods are necessary. In addition, 
viscous hypersonic and nonequilibrium flow problems are generally stiff and implicit methods are 
generally preferred over explicit methods. Some of the numerical issues for steady inviscid hypersonic 
blunt-body flow computations were addressed in our earlier paper [1]. A semi-implicit method and 
a fully implicit method for steady-state nonequilibrium flows were discussed in Yee and Shinn [2]. 
A basic study on numerical methods for unsteady inviscid nonequilibrium flows was presented in 
LeVeque and Yee [3]. The objective of this research is to efficiently extend and improve as well as to 
present an unified formulation of the existing implicit high-resolution shock-capturing schemes [4-6,2] 
for multidimensional compressible Euler and Navier-Stokes equations in the hypersonic, perfect and 
equilibrium real gas flow regimes. 

1 An abbreviated version will appear in the proceedings of the BAIL V Conference, June 20-24, 1988, Shanghai, China. 
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3 Research Scientist , 

4 Research Scientist, Theoretical Aerodynamics Division, this work was performed while on leave as an Ames Associate 

at NASA Ames Research Center, Moffett Field, CA 94035 USA. 


The improved schemes are based on a class of implicit Total Variation Diminishing (TVD) type 
algorithms originally designed for transonic and supersonic multidimensional Eider and Navier- Stokes 
equations [4-6]. The extended conservative shock- capturing schemes are spatially second- and third- 
order, and are and fully implicit. They can be first- or second-order accurate in time and are suitable 
for either steady or unsteady calculations. In addition, the current unified formulation allows the 
inclusion of the MUSCL-type approach [7] in conjunction with a local characteristic approach [ 24 , 6 ] 
or flux- vector splittings [8] (see section II for an explanation). For the present study, particular 
emphasis is placed on second-order implicit time-accurate high-resolution algorithms. The algorithms 
are formulated in finite volume and pseudo finite volume forms which, for certain physical problems 
and grid distributions, can enhance stability and convergence rate for highly clustered or skewed grids 
and require only a slight modification from the form originally presented in Yee and Harten [5] for 
generalized geometries. It is emphasized here that the use of the term TVD-type schemes pertains to 
the property of the algorithm as applied to one- dimensional nonlinear scalar hyperbolic conservation 
laws or constant coefficient hyperbolic systems in a semidiscrete sense. Theoretical justification of 
the proposed fully discretized schemes on the preservation of TVD property for the general nonlinear 
scalar hyperbolic conservation laws is under investigation. Moreover, the high-resolution property of 
these schemes for multidimensional nonlinear systems of hyperbolic conservation laws is evaluated 
by numerical experiments. In particular the following numerical issues are addressed: 

1. Some numerical aspects of TVD-type schemes that affect the convergence rate for hypersonic 
Mach numbers find real gas flows but have negligible effect on low Mach number or perfect gas flows 
are identified. 

2. The performance of the various linearized implicit forms of the proposed schemes similar to the 
transonic flow study [4-6] is reexamined for hypersonic flows. 

3. The behavior of the proposed schemes with various temporal differencing but similar spatial 
discretization for inviscid and viscous flows is investigated. Studies indicated that their behavior in 
terms of stability and convergence rate is quite different between viscous and inviscid flows. However, 
with the proper choice of the temporal discretization and suitable implicit linearization, these schemes 
are fairly efficient and accurate for very complex two-dimensional hypersonic inviscid and viscous 
shock interactions. 

4. The relative efficiency and accuracy of typical TVD-type schemes [9-11] for shock wave compu- 
tations are examined. A comparative study on steady and unsteady flows reveals that the proposed 
class of TVD-type schemes, in particular, for equilibrium real gas and nonequilibrium flows, pro- 
duces just as accurate shock resolution and yet requires less operations count than most other TVD 
schemes (e.g., higher-order Godunov [10], Osher & Chakravarthy [9], and TVD flux-vector splitting 
approaches [8]). 

In the following section, the generalization of Roe’s approximate Riemann solver and flux- vector 
splitting for real gases are reviewed. A description of the two-parameter family implicit shock- 
capturing scheme, the various enhancements on numerical stability and convergence rate for hy- 
personic flows, and the behavior of the scheme for inviscid and viscous flows are discussed in the 
subsequent sections. To illustrate the performance of the schemes for complex hypersonic flows, some 
representative numerical examples are also discussed. 

II. Description of the Numerical Algorithm 

The conservation laws for the two-dimensional Navier- Stokes equations can be written in the form 

dU dF dG _ J_ 
dt ^ dx ^ dy Re 


dF v dG v 
dx ^ dy 


(la) 


2 
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where U = [ p, m, n, e ] T , F = [ pu, mu + p, nu, eu + pu ] r , (? - [ pv, mv, nv + p, ev + pv ] , F v 
= [ 0, r«, r xy , / } T , and G v = [ 0, r xy , r yy , g } T . Here p is the density, u and v are the velocity 
components, m - pu and n - pv axe the x- and y-components of the momentum per unit volume, p 
is the pressure, e = p[e + |(u 2 + v 2 )] is the total energy per unit volume, and € is the specific internal 
energy. For a perfect gas, we also have 


r xx = u x - 2u y )/3, 

(ib) 

r xy = p{u y + v x ), 

(lc) 

T yy ~ M( — 2u x + 4v y )/3, 

(Id) 

. , da 2 

f = ur xx + VT xy + pPr (7 - 1) — , 

(le) 

g r= UT xy + V T yy + ^P r (l ~ 1) ’ 

(If) 


where for example, u x is defined as du/dx. The dynamic viscosity p is given by Sutherland’s formula. 
The Reynolds number is Re, the Prandtl number is Pr, the sound speed is a, and the ratio of specific 

heats is 7 . 

Under a generalized coordinate transformation, £ = £(*,y) and r) = T}(x, y), equation (1) can be 
written in a form which maintains the strong conservation-law form as 


du ap dG 1 \dF\, dGJ 
~dt + ~dl + dr) ~ K + d n r 


( 2 ) 


where V - UfJ , F = (£*F + £ y G)/J, G = {t] x F + rf y G)/J , F„ = {£ X F V + £ y G v )/J , G v = (Tx-Fv + 
r] y G v )/J , and J = f x rj y - £ v »}*, the Jacobian transformation. Let >1 = dF/dt/ and P = dG/dU. 

Then the Jacobians 4 = dF/dU and P - dG/dU can be written as 


A = (£ X A + ( y B) (3a) 

B = (i, x A + rj y B). (3b) 

In this study the thin-layer Navier-Stokes approximation is assumed by dropping all the d(-)/d£ 
derivatives in the viscous terms. Also, stability and convergence rate viscous results are for a perfect 
gas and laminar flows with adiabatic wall conditions. 


2*1. Riemann Solvers 

Here the usual approach of applying the one- dimensional scalar TVD schemes via the so called 
Riemann solvers for each direction in multidimensional nonlinear systems of hyperbolic conservat ion 
laws (see for example reference [5,12]) is used. This approach is best suited for orthogonal or nearly 
orthogonal grids. The eigenvalues and eigenvectors of the Jacobian matrices A and P are used in 
approximate Riemann solvers. Given two states whosedifference is A17, Roe [13] obtained an average 
A in the ^-direction, for example, satisfying A P = AAU for a perfect gas. The generalization by 
Vinokur [14] for an arbitrary gas involves the pressure derivatives \ = {dp/ dp)~ and k = {dpjde) p , 
where Z = pe. The relation c 2 = \ + /ch then gives the speed of sound, where h = c + p/p. Introducing 
H = h + (u 2 + t> 2 )/2, Vinokur found the same expressions for u,v and H as for the perfect gas, and 
that T and k must satisfy 
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\Ap + kAF= Ap. (4) 

Unique values of x and « are obtained by projecting the proper averages of the values for the two 
states into this relation (see references [14,11,12,15] for the exact formulas). 

Flux-vector splitting methods divide the flux F into several parts, each of which has a Jacobian 
matrix whose eigenvalues are all of one sign. The approach by Steger and Warming [16] made use 
of the relation F — AU , valid for a perfect gas. Van Leer [17] constructed a different splitting in 
which the eigenvalues of the split-flux Jacobians are continuous and one of them vanishes, leading 
to sharper capture of transonic shocks. Vinokur and Montagne [18] showed that the expressions for 
both these splittings can be generalized to an arbitrary gas by using the variable 7 = pc 2 jp, and 
adding to the split energy flux a term equal to the product of the split mass flux and the quantity 
e - c 2 /[7(7 - 1)] (see references [18,11,12,15] for the exact formulas). 

The current study on the shock resolution of the various schemes [1,4-6, 9] for two-dimensional 
steady-state blunt-body inviscid computations indicates similar trends as the one dimensional study 
[11], The main issue appears to be their relative efficiency. Due to extra evaluations per dimension 
in the curve fitting between the left and right states in a real gas for the van Leer formulation, 
additional computation is required for the van Leer type schemes than the Harten and Yee [4,5,19], 
and Yee [6] types of TVD schemes. Here van Leer type schemes refer to the use of the MUSCL 
approach [7] (see section 2.2 for an explanation) in conjunction with the Roe type approximate 
Riemann solver [13] or flux- vector splittings [8] (hereafter the latter methods are referred to as the 
TVD flux-vector splitting methods). Moreover, for steady-state applications, implicit methods are 
preferred over explicit methods because of the faster convergence rate. In addition, it is easier to 
obtain a noniterative linearized implicit operator for the Harten and Yee, and Yee type schemes 
than for the van Leer type schemes. Furthermore, unlike flux-vector splitting approaches, implicit 
methods employing the Roe type approximate Riemann solver (non-MUSCL or MUSCL) with first- 
order implicit operators do not require the Jacobian of the F ± and G* fluxes. Here F ± is the 
portion of the flux with positive/negative eigenvalues. In many instances, the Jacobians of these 
fluxes are relatively difficult or expensive to obtain, in particular for nonequilibrium flows. A similar 
difficulty applies to the MUSCL formulations via the Roe type approximate Riemann solver if a 
spatially second-order implicit operator is desired. For these reasons, the linearized implicit versions 
of Harten and Yee [4] and Yee [6] are preferred over the van Leer type schemes. Consequently, 
numerical studies on the extension of the former schemes to hypersonic flows are emphasized. Some 
of these points will become more apparent when an unified formulation of these implicit methods is 
presented in section 2.2. An unified formulation of the corresponding explicit schemes can be found 
in reference [12]. 

2.2 Description of the Implicit TVD schemes 

In the application of TVD-type schemes for viscous flows, the physical problems considered here 
are assumed to be inviscid dominated in the sense that moderate or strong shock waves are present 
in the flow field such that high- resolution shock- capturing techniques are required. Thus the numer- 
ical procedures used here for the compressible Navier- Stokes calculations are a second-order central 
difference approximation for the diffusion terms and TVD-type schemes for the inviscid part of the 
Navier-Stokes equations. The question of whether the present numerical dissipation term (due to 
the TVD-type terms) has an adverse effect on the true viscosity terms in the boundary layer region 
is not known at this point. What we can conclude from the current study is that the portions of the 
solution slightly or far away from the boundary layer are quite accurately simulated. 

The two-parameter family of explicit and implicit high-resolution schemes presented here is based 
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on a semi-discrete methodology and on the one-parameter family of TVD-type algorithms developed 
in references [19,4-6]. The idea is to use the same spatial discretization as references [19,4-6] for 
the spatial derivatives and to use the two-parameter family of linear multistep methods for the time 
derivatives. The original one-parameter family of TVD-type schemes is a subset of the two-parameter 
family of algorithms. Mathematical analysis similar to that in [19,4-6] for the current larger family 
of schemes is under investigation. For a particular chosen time differencing, these schemes are 
TVD for the one- dimensional constant coefficient hyperbolic equations. Also the MUSCL approach 
in conjunction with the Roe-type approximate Rieinann solver [2] and TVD flux-vector splitting 
methods [8] falls nicely into the present frame work. In other words, the present formulation includes 
a larger class of spatial as well as temporal discretization than in references [19,4-6]. 

Also, the formulation is in finite volume and pseudo finite volume forms which can enhance stability 
and convergence rate for highly clustered or skewed grids and is a slight modification from the form 
originally presented in Yee and Harten [5] for generalized geometries. For fairly uniform or mildly 
clustered grids, the present finite volume and pseudo finite volume forms behave the same as in 
reference [5] for inviscid flows. This is in contrast to the study of Takakura et al. [20] which claimed 
that their modified form [20] is the correct finite difference formulation for generalized geometries. 
A comparison between Takakura et al. [20] and reference [5] on the same fairly uniform curvilinear 
grid for a blunt-body computation shows no noticeable difference in resolution. 

Without loss of generality, the two-parameter family of implicit schemes for the Euler equations 
(F„ = G v = 0) is presented here. For general Navier-Stokes equations, the appropriate three-point 
central differences of the viscous Jacobian terms should be added to the implicit operator and a 
central difference approximation for the diffusion terms should be added to the explicit operator. 

Let At be the time step and let the grid spacing be denoted by Af and A rj such that £ = ;'A£ 
and T] = fcAq. Also let A« = and A" = then a two-parameter family of explicit and implicit 
TVD-type algorithms in generalized coordinates for two-dimensional systems (1) with F„ = G v = 0 
can be written as 


+ 


— fF n+1 

1 + u; 1 


?n+ 1 


,1 + 1 


d-^fSn 

l F j+b,k 


1 + W 




(l-*) A" 


1 + u> 


[Gj. 


fc +1 




+ 




a u"; 1 . 

l + w J ' k 


(5) 


Here A 17" = U 


n+ 1 


77 n 

U jf 


The functions F j+ 1 k and G j k+ i are the numerical fluxes in the £- 


j t k — u j,k ■ - 

and 77-directions evaluated at ( j + |,Ar) and ( j, k 4- 5), respectively. This two-parameter family 
of algorithms contains first- and second- order implicit as well as explicit schemes. The scheme is 
temporally second-order if 0 = u; + \ and first-order otherwise. When <9/0, algorithm (5) is an 
implicit scheme. In this paper, only the temporally first-order backward Euler (6 = 1, w = 0) and the 
temporally second-order three-point backward differentiation (8 — 1, u; = 1/2) time-differencing are 
investigated. Detailed formulation and numerical studies for algorithm (5) with = 0 for transonic, 
supersonic and hypersonic flows can be found in references [4-6,12,1,11,21,22]. The current study 
shows that, for viscous steady and unsteady flows, the temporally second-order implicit algorithm 
(8 - 1, u? = 1/2) appears to be slightly more stable and efficient than the temporally first-order 
implicit algorithm (8 — 1, = 0). 


The spatial accuracy of equation (5) depends on the form of the numerical flux functions. There 
exists many ways to achieve higher-order spatial accuracy and at the same time have TVD-type 
properties. Here two of the ways are discussed. The first is due to Harten [19], Roe [23], and 
Yee-Roe-Davis [6,24,25], and the second, sometimes referred to as the MUSCL approach, is due to 
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van Leer [7]. Following the same nomenclature as in references [11,12], hereafter, we refer to the 
first way as the non-MUSCL approach. Besides the different temporally implicit discretization, the 
combination of the two Riemann solvers discussed in section 2.1 and higher-order spatial differencing 
considered yields three different types of spatial differencing for the nonlinear system (1): namely, 
non-MUSCL and MUSCL approaches using an approximate Riemann solver, and a MUSCL approach 
using flux- vector splittings (TVD flux-vector splitting methods). 

Non-MUSCL Approach Using an Approximate Riemann Solver. The numerical flux funct ion F ,+ M 
for a non-MUSCL type approach, together with the local characteristic approach [4-6] (Roe type of 
approximate Riemann solver) in a finite volume formulation, can be expressed as 


Fj + 1 ,k — 2 ( j") . i ( Fj,k + + ( j) + i^'* + + Rj+i$ J+ i j Jj+ 5 


(6a) 


The corresponding pseudo finite volume formulation will be discussed in section 2.4. Here Rj + 1 

is the eigenvector matrix for dF/dU evaluated at some symmetric average of Uj,k and Uj+ 1 ,* (for 
example, Roe average [13] for a perfect gas and generalized Roe average of Vinokur [14] for real 
gases). The values 




1 _ 1 / 1 1 \ 
•7j+ 1 2\Jj+l,k jj,k ) 


(6b) 


Also (h)j+i - 


are defined, for example, as 


VVW + (¥j 


->'+a 


and {h) j+ i 




+ (tV 


i i+\ 








used in R , i 

J-r 2 


(6c] 


The values £ x , £,,, tj x and rj y are evaluated by three-point central differences. Similarly, one can 
define the numerical flux G } in this manner. 

Here the form of 4-1 can be divided into two types: (a) a spatially second-order symmetric 
TVD-type scheme [6,24,25] in which the numerical dissipation terms are independent of the sign of 
the characteristic speeds and (b) a spatially second-order upwind TVD-type scheme [19,5] in which 
the numerical dissipation terms depend on the sign of the characteristic speeds. 

s 

The elements of the $ l4 i in the ^-direction denoted by (<f>‘ j+ i) for a spatially second-order 


symmetric TVD-type scheme [6,12] are 




(7a) 


The value a 1 , i is the characteristic speed a 1 for dF/dU evaluated at the same symmetric average 

J+ 2 

between U it k and Uj+i >fe . The function ip is 


xp(z) 


1*1 

(z 2 + i 2 )/2^i 


1*1 > *i 
1*1 < 


(7b) 


Here xp(z) in equation (7b) is an entropy correction to |z| where 6 x is a small positive parameter. 
Por problems containing simple unsteady shocks, S\ is set to zero in most of the computations 
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since entropy- violating phenomena occur only for steady or nearly steady shocks. For steady-state 
problems containing strong shock waves, a proper control of the size of 5, is very important, especially 
for hypersonic blunt-body flows. The choice of 6 X is also highly dependent on the Mach number and 
geometry of the physical problem. For the current numerical examples, the parameter S x is set to be 
a function of and aj,. See reference [12] or section III for a discussion. 

Examples of the ‘limiter’ function Q‘ j+ i can be expressed as 


<?^ + i = minmod(</_i,a' + i) + minmod^+i,^*) - + i 

Q l j+ 1 =minmod(a'._i,af + i,a'. + j) 

Q l j+ i =minmod[2a'_i,2a^ + i,2a'. + i,-(a^_i +«' i+ i)]- 


(7c) 

(7d) 

(7e) 


The minmod function of a list of arguments is equal to the smallest number in absolute value if the 
list of arguments is of the same sign, or is equal to zero if any arguments are of opposite sign. Here 


a 1 , i are elements of 
■?+ 2 


Q J+ i = R - Uj'k)- 


(8) 


The elements of the in the ^-direction denoted by (^ + i) f° r a spatially second-order 

upwind TVD-type scheme [19,5,12] are 



(9a) 


where 


7* + i = 







Examples of limiter function g\ used in calculations are 


(9b) 


g'j = minmod( a* _ i , i ) 

»; = ("'+ i"!!-! + / ("i+i + ">-s) 

= {«' - J [<«' ♦ 5 ) ! + «] + S H- { ) | ! + «] } / K 1 : >’ + K- i : > ! + 24 
9 l j = minmod(2a'_i,2aJ + i,haJ_i + a^i)) 


g l j ; = 5 • max 


0,min(2]a' j+ i|,S • a*_i),min(|a* + i|,2S • «j_i) 


; S = sgn(a^ + i). 


(9c) 

(9d) 

(9e) 

(9f) 

(9g) 


Here 6 is a small parameter. In practical calculations 10 7 < 8 < 10 5 is a commonly used range. 


MUSCL Approach Using an Approximate Riemann S olver . The numerical flux function F J+ i for 
a MUSCL type approach, together with Roe type of approximate Riemann solver, for an upwind 
scheme as described in Yee [26] and Yee & Shinn [2] can be expressed as 
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The values F(Uf + ^) and F{Uf + ^) are the flux function F evaluated at Uj * j and Uj^ * respectively, 
with 


P* j = U,+1* - \ [(1 - + (1 + 3) A j+i] 

(10b) 

v&j = Vij. + i[(l - flATTj + (1 + 

(10c) 


where rj discussed below, is a parameter to control the spatial accuracy of the scheme. The limiters 
and can be expressed as 

J-r 3 J~r 3 

A J+ i = minmod(A^ + x,^A J _i) ( 10d) 

A j+ ^ = minmod(A j+ i,/3A j+ i) (lOe) 

with , 

minmod(.r,/3j/) = sgn(x) • maxi 0,min[|a:|,sgn(a:) • (3y] ^ (lOf) 

and (3 = y5|, where A J+ i = U j+ i >fc - U Jik - For fj - -1, A j+ i and A J+ i can be the same limiter as 
(9) except the arguments Me now the A J± i instead of a J± i. 

The vector j? J+ i is the eigenvector of A evaluated at some symmetric average of an d 1 ; 
i.e., 




j+i 1 


(lOg) 


and the elements of are 


= ^( 2 ' + 1 ) 5 ' + 1 


j+: 


j+ : 


where again a^ +i and af +1 are evaluated at some symmetric average of and Uf +k and 


j+ 




(lOh) 


(lOi) 






J+ 


(ioj) 


Here the spatial order of accuracy pertaining to the scheme with the limiter not present (i.e., 
remove the tildes from equations (10b,c) ) is determined by the value of rj. For example, if rj = -1 
the scheme is fully upwind and if rj = 0 it is Fromm’s scheme. When rj = 1/3 the scheme is third-order 
and when rj = 1 it is the regular three-point central difference scheme. 

MUSCL Approach Using Flux-Vector Splittings: The numerical flux function F j+ 1 k for a MUSCL- 
type approach, together with flux-vector splittings [8] referred as the TVD flux-vector splitting 
method in this paper can be expressed as 
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! 


w = i> + a+(u ‘ 


L ' 

j'+3 ; 


(11) 


where f’ ± (t/ R ’^) are evaluated using either the Steger- Wanning type [16,18] or van Leer type [17,18] 
flux- vector splittings. The vectors U* + * and Uf + * are the same as in equations (10b, c). The quantity 
F~(U R x) is the portion of the flux F with negative eigenvalues evaluated at Uf + k . 

The operations count between (6)-(9) and (10,11) is within 30% for a perfect gas. However, 
due to an extra evaluation per dimension in the curve fitting between the left and right states m 
an equilibrium real gas for (10,11), additional computation is required for the MUSCL approach. 
The slight advantage of (10,11) over (6)-(9) is that (10,11) can be spatially third-order accurate. 
However, experiences with the third-order case (rj = 1/3) do not show a very visible improvement 
over the second-order case for problems with discontinuities. Part of the reason is that all TVD-type 
schemes reduce to first-order at points of extrema regardless of the order of accuracy at smooth 
regions Also, because of the similarity in shock resolution and yet higher operations count for 
real gases and nonequilibrium flows of the MUSCL over the non-MUSCL approach using Roe type 
approximate Riemann solver, efforts are concentrated only on the non-MUSCL formulation. At 
present no outstanding advantages or disadvantages between these formulations for a perfect gas 
have been observed. Further investigation is required along this line before a clearer comparison can 

be drawn. 


2.3 A Conservative Linearized Implicit Form For Unsteady and Steady-State Calcula- 
tions 

To solve for U n+l in (5) one normally needs to solve a set of nonlinear algebraic equations itera- 
- tively. One way to avoid this is to linearize the implicit operator and solve the linearized form by 

other means. Following the same procedure as in references [4-6], a conservative linearized alter- 
I nating direction implicit (ADI) form of (5) for the numerical fluxes (6) and (10a) can be written 

as 


1 + U) J+3 ‘ k 1 + W -* 3 ’ 


where 


E* = - 

+ 


A* 

1 + 

U) 




A" 

1 + u> 


a u: 


n- 1 




1 + u> J ,fe ’ 

A"0 


Hr 


1 -f OJ 2 1 + U> i' k 2 \ 

E n = A f/ n ; U n+l = U n + E n , 




£" = E\ 


The nonstandard notation 


1 r r 

2 


Hj + i k E* - -[Aj + i' k Ej +l k - D^ + i fc F*] 


p-in /~rn 

(12a) 

(12b) 

(12c) 

(12d) 

(12e) 

(12f) 
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is used, and ft* 

7 i 


j+\,k 


E * , n7. ji-E can be taken as 


n j + i ifc £ * = 5j+i,fcdiagN 



(12g) 

(12h) 


Here -Aj+i,*, Bj t k + i are Jacobians of F and G evaluated at (j + 1 , fc) and (j, fc + 1). The value 
E” k = C/j 1 * 1 - C/J 1 ^ . The expression diag(r / ) denotes a diagonal matrix with diagonal elements z l . 

The nonconserverative linearized implicit form suitable for steady-state calculations [4] is also 
considered- Numerical study indicated that the latter form appears to be slightly less efficient in terms 
of convergence rate than the linearized conservative form (12). This conservative linearized implicit 
operator as well as the nonconservative linearized implicit operator (both suggested in reference [4]) 
for 6 = 1, a? = 0 was renamed two year later by Barth [27] as the approximate Jacobian linearization. 
To compute (12g,h), a triple matrix multiplication of dimension (4 x 4) has to be performed at every 
grid point. For steady-state applications, one can simplify (12g,h) as 


tt] + i tk E = H E j+i,k ~ Ej' k ) (13a) 

^ T1 j: k+i E = M n I{ E j,k + 1 - (13b) 

The scalar values and M ^ are 

~ max ) (13c) 

Mrj - max0(a^ + i) (13d) 

and J is the identity matrix. Note that (13a,b) involve scalar multiplication only. The solution using 
(13) is still second-order (or third-order) accurate after it reaches steady-state. Other linearized 
implicit forms can be found in references [4-6]. 

2.4 General Assumptions and Limitations on the Numerical Studies 

The present study is by no means an exhaustive investigation. There are additional elements and 
parameters (other than the ones considered here) in the algorithm itself as well as in the physical 
problem, such as flow type and geometric complexity, that can affect or interfere with the performance 
of the numerical scheme. Even within the numerical issues listed in section I, the study is limited 
to a sampling of the parameter range for the time-step size or OFL number and the form of in 
(7b). In particular, various strategies to speed up and stabilized the start-up solution from freestream 
conditions for steady computations have not been investigated. What is discussed here is intended to 
give interested readers some guideline for the use of the algorithm and to shed some light for further 
study and improvement of the scheme and the development of better ones. All of the numerical 
studies discussed in the subsequent sections rely on the following assumptions and considerations: 

1. Although the recommended finite volume formulation (6) closely mimics the regular finite 
volume formulation for two dimensions, the results obtained in this report used a slightly different 
formulations than (6). In particular, three formulations (hereafter referred to as the pseudo finite 
volume formulations) for the non-MUSCL schemes were investigated and are as follows 
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(14a) 


F j+ I ,k 2 




(€.)i+|( J? M + ^J+i. fe ) + Uv)j+i(^J. fe + ^i+ l > fe ) + 

with the corresponding quantities (( x ) j+ ii J j+ \ and (*i)j+£ of equations (6b, c) express as 

(£*) i+ i = ^[U-)j+i.fc + = + (Ub) 

(14c) 


(^1 )j+4 - 


(&W 


yUxj^i + (Oj+j 


and 


FV 


j+£,fc “ 


(^ + ^) jifc +(^ + ^G) i+1> 


'i+i 


(15) 


and 


^j+ 1 .fc 


•^j+i.fc + ^J, fc + 


A +s . 


(16) 


Here J + i and ifci in equations (15) and (16) are the same as (14b, c). For highly skewed grids and 
nonuniform flows, equations (6) and (14) are preferred over (15) and (16). However, (14) and (15) 
do not preserve freestream whereas equations (6) and (16) do. All of the results present in section V 
use (15). One of the blunt-body cases was rerun with equation (6) and (14)-(16) and no noticeable 
difference was observed. We expect all of the conclusions on the behavior of (14)-(16) to be carried 
over to equation (6), since all of the examples use mildly clustered yet quite regular and nearly 
orthogonal grids. 

In two dimensions the present pseudo finite volume formulations can be made ‘truly finite volume 
by a slight modification of equations (14)-(16); i.e., on the treatment of k { and J j+ i. However, 
the situation is different in three dimensions where finite volume formulations depart from finite 
difference formulations. See reference [28] for a discussion. 


2. The numerical results and conclusions are for the non-MUSCL approach and for the particular 
flow type and geometry specified with a sampling of a narrow range of Mach numbers and time steps. 
Results for viscous flow calculations are based on the shock wave dominated thin-layer Navier-Stokes 
equations for lanunar flows. 

3 The numerical boundary condition treatments are the same as in references [5,29,30] for the 
inviscid flow and as in reference [21] for the viscous flows. Studies [31] showed that proper treatment 
of numerical boundary conditions has a major impact on the stability and convergence rate of the 
scheme. Therefore the types of boundary condition treatment used here reflect on the performance 
of the stability, accuracy and convergence rate of the present algorithm. 


4. For steady-state computations, the convergence rate not only depends on all of the elements and 
parameters (to be discussed shortly), but more importantly also on the type of grid associated with 
the computation. Studies show that, in general, a coarse nearly uniform orthogonal grid converges 
1-3 times faster than a similar finer grids, and possibly an order of magnitude or more faster than 
highly clustered or skewed grids. What will be presented in section V represents fairly uniform to 
mildly clustered grids. Most of the grids used for the numerical study were not very coarse; thus the 
number of iterations quoted is naturally higher than its coarse grid counterpart. 
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5. For the non-interfering blunt-body flows, the convergence rate and behavior of the symmetric 
and upwind TVD-type schemes are very similar. However, for the interfering blunt-body flows 
containing slip or shear surfaces, the upwind scheme produces sharper weak solutions. Consequently, 
all of the illustrations and conclusions discussed in this paper are for the upwind scheme using limiter 
(9c). Other limiters can produce sharper discontinuities but are not as robust as limiter (9c). See 
section V or reference 12 for a discussion. 

6. Research on real gas effects on the performance of the proposed scheme is only in the preliminary 
stages. All of the illustrations and conclusions for real gases are for invsicid non- interfering blunt- 
body flows. Study of viscous real gets flows is in progress. 

7. For steady-state computations using the backward Euler time differencing (0 = 1, u = 0), a 

local time stepping procedure similar to [30,32] was used. However, in comparing the convergence 
rate with the three-point backward differentiation time differencing (9 — 1, = 1/2) for the viscous 

flows, a constant time step was used. 

$ 

Other considerations such as reducing ADI factoring error, using multigrid, relaxation or conjugate 
gradient methods as an alternative to ADI, using local grid refinement to enhance resolution, etc., 
are also sources of improvement to algorithm (6-16). These items and the development of better 
algorithms are the subject of on-going research. 


III. Enhancement of Stability and Convergence Rate for Hypersonic Flows 

In reference [1], some elements and parameters which can affect the stability and convergence rate 
in high Mach number cases but have negligible effect in low Mach number cases for steady-state 
inviscid blunt-body flows were identified. The current study indicated that the same elements and 
parameters can affect the stability and convergence rate at hypersonic speeds for viscous computa- 
tions as well. They are as follows: ( 1 ) the choice of the entropy correction parameter , (2) the choice 
of the dependent variables on which the limiters are applied, and (3) the prevention of unphysical 
solutions during the initial transient stage. 

1. For Mach number ranging from 1.2 to 15, numerical experiments for one- and higher- 
dimensional unsteady flows containing unsteady shocks show that the second-order explicit TVD 
schemes [29,12,11] are insensitive to the entropy correction for 0 < 6\ < 0.1. In most cases 8\ = 0 
was used. For 0.1 < < 0.25, there is a possibility of improving stability in the sense of allowing a 

higher CFL number at the expense of a slight smearing of the discontinuities. However, for unsteady 
complex shock wave interactions, a small positive <5i or a variable 6\ (to be discussed) cam help 
stabilize the time-accurate implicit algorithm (12). 

For subsonic to low supersonic steady-state NACA 0012 airfoil computations [5,6], the resolution 
of the shock waves was found to be quite insensitive to 0.1 < 6\ < 0.125 and a constant value seems 
to be sufficient. However, for hypersonic flows, especially for blunt-body flows, a constant 8\ or 
a variable 5\ suggested by Harten and Hyman [33] was found to be insufficient, but a variable 8\ 
depending on the spectral radius of the Jacobian matrices of the fluxes is very helpful in terms of 
stability and convergence rate. In fact, a proper choice of the entropy parameter 8\ for higher Mach 
number flows not only helps in preventing nonphysical solutions but can act, in some sense, as a 
control of the convergence rate and of the sharpness of shocks and slip surfaces (or shear layer in 
viscous flows). The smaller the 8\ that is used, the slower is the convergence rate. The larger the 
tfi that is being used, the larger is the numerical dissipation being added. However, 8 j cannot be 
arbitrarily large. 
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For the present blunt-body steady-state calculations with Mach numbers M > 4, the initial flow 
conditions at the wall are obtained using the known wall temperature in conjunction with pressures 
computed from a modified Newtonian expression [34]. Also, for implicit methods a slow startup 
procedure from initial conditions [30] is necessary. Most importantly, experience indicates that if one 
sets <§i in equation (7b) as a function of the velocity and sound speed, i.e., 

(«i) J+ j = <(i« j+ ji + i» i+ ii + c j + i> (17s) 

=MK+}I + K+{I + '*+!> (17b) 

with 0.05 < 6 < 0.25, then blunt-body flows for 4 < M < 25 appear to be stabilized and nonphysical 
solutions Me less likely to occur. Equation (17) is written in Cartesian coordinates. In the case of 
generalized coordinates, the u and v should be replaced by the contravariant velocity components, 
and one half of the sound speed would be from the ^-direction and the other half would be from 
the q-direction. For implicit methods, it is very important to use (17) in V>(~) on both the implicit 
and explicit operators since in a two-dimensional hypersonic flow field consisting of a mixture of 
subsonic and supersonic regions with Mach number ranging from 0 to hypersonic speeds, an entropy 
parameter like (17) is nonzero in all of the regions. The entropy parameter (17) seems to work 
well for blunt-body flows but whether this is also the right choice for configurations other than a 
blunt-body shape is an open question* 

For unsteady hypersonic blunt body complex shock wave interactions, the entropy parameter 
(17) can help stabilize the time-accurate implicit algorithm. For most of the viscous and inviscid 
calculations shown, unless otherwise indicated, 6 is set to 0.125. 

2. Higher-order TVD schemes in general involve limiter functions. However, there are options 
in choosing the types of dependent variables when applying limiters for systems of hyperbolic con- 
servation law, in particular for systems in generalized coordinates. The choice of the dependent 
variables on which limiters are applied can affect the convergence process. In particular, due to 
the nonuniqueness of the eigenvectors R jJr i , the choice of the characteristic variables on which the 
limiters are applied play an important role in the convergence rate as the Mach number increases. 
For moderate Mach numbers, the different choices of the eigenvectors have a negligible effect on the 
convergence rate. However, for large Mach number cases, the magnitudes of all the variables at the 
jump of the bow shock are not the same. In general, the jumps are much larger for the pressures 
than for the densities or total energy. Studies indicated that employing the form R j+ i such that 
the variation of the a are of the same order of magnitude as the pressure would be a good choice for 
hypersonic flows. The form similar to the one used by Gnoffo [35] or Roe and Pike [36] can improve 
the convergence rate over the ones used in references [37,38]. In all of the computations shown the 
form Rj + 1 used is the same as in references [37,38] except for an extra factor of l/c^ + i • With this 
extra factor the variation of the a are in fact proportional to the pressure. Other forms of R J+ i 
have not been investigated. 

3. Due to the large gradients and to the fact that the initial conditions are far from the steady- 
state physical solution, the path used by the, implicit method can go through states with negative 
pressures if a large time step is employed. A convenient way to overcome this difficulty is to fix a 
minimum non-negative allowed value for the density and the pressure. With this safety check, the 
scheme allows a much larger time step and converges several times faster. 

In addition, since the Roe’s average state allows the square of the average sound speed c*. + i to lie 
outside the interval between c j and cj +1 , the average state c j + i might be negative even though Cj 
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and c* +1 are positive during the transient stage when the initial conditions are far from the steady- 
state physical solution. In this case, we replace c? +i by max(c* + i,min(e],c* +l )). This latter safety 
check is in particular helpful for the symmetric TVD algorithm (7). 


IV. Behavior of the Algorithm with Different Temporal Differencing 

It is emphasized here that since the method (12) is written in the ‘delta’ formulation, either 
the backward Euler (first-order) or the three-point backward differentiation (second-order) time 
discretizations require the same amount of storage and a similar operation count. Therefore, the main 
consideration between the two time- differencing methods is their relative stability and convergence 
rate. 

Inviscid Unsteady Flows : For inviscid unsteady flows, the explicit TVD-type methods [29,12,11] are 
more efficient than the second-order implicit method (12). Unless the inviscid problem is stiff, there 
is no advantage of employing an implicit method for inviscid unsteady flows. 

Inviscid Steady Flows : The backward Euler implicit method has a better stability and convergence 
rate than the three-point backward differentiation implicit method. Also a local time-stepping pro- 
cedure can speed up the convergence rate for the former time- differencing method whereas the same 
procedure has little effect on the convergence rate when compared with a fixed time step procedure 
for the latter time-differencing method. 

Viscous Unsteady Flows: Computations on the unsteady viscous flows mainly use the second-order 
time differencing since a larger time step can be used compared with the temporally first-order im- 
plicit method. Due to the highly clustered viscous grid used in contrast to their inviscid counterpart, 
solving a viscous unsteady complex shock interaction using an explicit TVD-type method is not 
practical due to its inherent time step restriction. In certain cases, the time step might be an or 
der of magnitude smaller than the implicit counterpart. A more detailed study of unsteady viscous 
hypersonic blunt-body flows with an impinging shock is reported in reference [22]. 

Viscous Steady Flows : At present there is no detailed viscous steady flow study comparing the first- 
order time differencing using a local time-stepping approach with the second-order time differencing 
using a constant time step approach. But the general trend is that the second-order time differencing 
has slightly better stability and convergence rate than the former one. In particular, a summary using 
a fixed time step approach comparing the two time- differencing algorithms is discussed in section V. 


V. Numerical Results 

The various numerical aspects discussed in sections III-IV are complimented by a variety of steady 
and unsteady, viscous and inviscid hypersonic blunt-body flow computations in this section. Six 
types of blunt-body test cases are illustrated in figures 1-11. Test cases 1 and 2 are inviscid, perfect 
and real gas, non-interfering blunt-body flows. Test case 3 is a steady inviscid, perfect gas blunt-body 
flow with an impinging shock. Test cases 4-6 are viscous steady and unsteady perfect gas blunt-body 
flows with and without impinging shocks. 

Comparison Amo ng the Various Linearized Implicit Methods : For the implicit operator, numerical 
experiments show that the linearized conservative form (12) converges slightly faster than the lin- 
earized nonconservative form [4] for both viscous and inviscid flows. It seems also that when the 
freestream Mach number increases, the convergence rate of the linerarized conservative form (12) 
is better than a simplified version which replaces and of (12g,h) by max ( ip{a l j+ i) 
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and max; xp{a l . x ) times the identity matrix (equation (13)). This is especially true for viscous flow 
computations. Due to the experience gained by the transonic and the inviscid hypersonic study, no 
detailed computations using the linearized nonconservative form were performed for viscous steady 
flow. All of the results and discussions for the viscous computations are based on the conservative 

linearized form. 

Another area of investigation is that for viscous computations, the Jacobian of the viscous terms 
on the implicit operators are rather expensive to compute. To maintain the spatial order of accuracy, 
for sure these terms are needed for unsteady flows. Whether the omission of these terms has a major 
impact on the stability and convergence rate of the algorithm for steady-state calculations is not 
known. Therefore, an investigation has been made on the difference in the convergence rate for the 
algorithm with or without the viscous terms in the implicit operator. A brief summary is me u e 
in the following subsection. 


Choice of Limiters: Unlike flows with transonic and low supersonic shock waves, problems containing 
strong hypersonic shock waves are more sensitive to the treatment of limiters. I sing the more 
diffusive limiter (7c) or (9c) turns out to be more stable than other more compressive limiters. In 
terms of shock resolution for both the symmetric and upwind TVD-type of schemes the sequences 
written in equations (7c)-(7e) and (9c)-(9g) are in order of increasing accuracy On the other hand, 
these sequences are in order of decreasing stability and convergence rate. The more compressive 
limiters like (9f) and (9g) have a very low stability and slow convergence rate for steady flows. e 
same conclusion applies for unsteady flows where the more compressive limiters have a very restricted 
time step limit. From our experiences, it is not advisable to use (9f) and (9g) for complex steady 
shock wave interactions. In particular, limiter (9g) should be used only for the linear fields (i.e., or 
the u and v characteristic fields in the x- and y-direction respectively). See reference [12] for more 

details. 

Convergence Rate of Explicit and Implicit TVD-type Schemes for Rea l Gas Flows: The five differ- 
ent second-order TVD methods previously studied [11] in one dimension yield very similar shock 
resolution for the blunt-body (non-interfering case) problem. In particular, for an inviscid blunt - 
body flow in the hypersonic equilibrium real gas range, the explicit second-order Harten and Yee, 
and Yee-Roe-Davis type TVD schemes [6,24,25] using the generalized approximate Riemann solver 
[13] produce similar shock-resolution but converge slightly faster than an explicit second-order van 
Leer type scheme using the generalized van Leer flux- vector splitting [11]. 


The freestream conditions for the current study are - 15 and 25, p » - 1.22 x 10" 3 NJm , 
p = 1 88 x 10' 2 kg/m 3 , and T«, = 226°K. Figure 1 shows half of the 61 x 33 grid used for the blunt- 
body problem. For the = 25 case, the shock stand off distance is at approximately fourteen 
points from the wall on the symmetry axis. The relaxation procedure for the explicit methods 
employs a second-order Runge-Kutta time discretization with a CFL of 0.5 (solution not shown). 
The parameter 6 is set to a constant value of 0.15. Pressure and Mach number contours converge 
and stabilize after 3000-4000 steps but the convergence rate is much slower for the density (with a 
2-3 order of magnitude drop in I 2 -norm residual). The bow shock is captured in two to three grid 
points. The curve fits of Srinivasan et al. [39] are used to generate the thermodynamic properties of 

the gas. 


The same flow condition was tested on the implicit scheme (12) and the convergence rate was 
found to be many times faster. Figures 2 and 3 show the Mach number, density, pressure and k 
contours computed by the linearized conservative ADI form of the upwind scheme (12) with the 
first-order backward Euler (6 = 1 and u, = 0) for Mach numbers 15 and 25. Figure 4 shows the 
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slight advantage of the convergence rate of the linearized conservative implicit scheme (12) over the 
linearized nonconservative implicit scheme (with 0 = 1, u; = 0 find a fixed CFL of 15) suggested in 
reference [4], The convergence rate and shock resolution for the symmetric TVD-type scheme (12) 
behave similarly. For — 15 case, the L 2 -noTm residual stagnated after a drop of four orders of 
magnitude. 

In general, for a perfect gas with 10 < < 25 and a not highly clustered grid, steady-state 

solutions can be reached in 600-800 steps with 12 orders of magnitude drop in the Z^-norm residual. 
However, the convergence rate is many times slower for the real gas counterpart. Figure 5 shows the 
convergence rate for a perfect gas compared with a real gas computation with a fixed CFL of 50. Note 
that the scale of the ordinates used in figure 5 for the perfect gas and the real gas are not the same. 
The freestream conditions for the real gas study are the same as figure 3. An important observation 
for the behavior of the convergence rate for the Mach 15 real gas case is that the discontinuities of 
the thermodynamic derivatives which exist in the curve fits of Srinivasan et al. [39] might be the 
major contributing factor. This is evident from figures 2d and 3d and from a comparison with the 
convergence rate for the perfect gas. Another contributing factor is that the curve fits are accurate 
only for temperatures up to 6000°K. Since the temperature in this case is slightly above 6000°K, 
there is an uncertainty in the accuracy of the computed results. Further improvement of the existing 
curve fitting procedure is needed. 

Inviscid Impinging Shock Computations : Figures 6 and 7 show the schematic of the computational 
domain, the Mach contours and Z^-norm residual computed by the implicit upwind scheme (12) (with 
9 = 1, a; = 0) of an inviscid shock- on- shock interaction on a blunt body with radius Ri and thickness 
D — 2 Ri in the low hypersonic range. Higher inviscid hypersonic Mach number computations using 
the proposed scheme are also possible but are not shown here. Some viscous and inviscid studies on 
flow fields of this type were reported in references [34,40-42]. This flow field is typical of what may be 
experienced by the inlet cowl of a hypersonic aerodynamic vehicle. The freestream conditions for this 
flow field are the freestream Mach number M = 4.6, the freestream temperature T <*> — 167°K, and 
7 = 1.4 for a perfect gas. An oblique shock with an angle of 20.9° relative to the free stream impinges 
on the bow shock. Various types of interactions occur depending on where the impingement point is 
located on the bow shock. As shown by the Mach contours ranging from 0 to 4.55 in increments of 
0.05, the impinging shock has caused the stagnation point to move away from its undisturbed location 
at the symmetry line. The surface pressures at the new stagnation point can be several times larger 
than those at the undisturbed location of the stagnation point. In addition, a slip surface emanates 
from the bow shock and impinging shock intersection point and is intercepted by a shock wave which 
starts at the upper kink of the bow shock. The interacting shock waves and slip surfaces are confined 
to a very small region and must be captured accurately by the numerical scheme if the proper surface 
pressures are to be predicted correctly. The 77 X 77 grid used and the convergence rate computed 
by the implicit scheme (12) are shown in figure 7. Though the pattern of the flow is significantly 
more complicated than for the previous cases, the convergence rate remains quite satisfactory. As 
shown in figure 6 at the inflow, all of the inviscid and viscous interfering blunt-body computations 
start with the appropriate freestream and oblique shock wave conditions as boundary conditions. 

Viscous Steady Computations With or Without Impinging Shock: To keep the study tractable only 
two types of physical flow fields were chosen. The first is the viscous hypersonic blunt-body flow 
at A/qc = 8.03 and T <*> — 122. 1°K with a laminar Reynolds number of 387,750 based on the body 
diameter for ideal gases. The second problem (with the same flow conditions) is similar to the 
inviscid shock-on-shock interaction where an oblique shock impinges on the bow shock of the blunt 
body at an angle of 19° relative to the free stream. A more detailed flow field computation on the six 
types of shock patterns categorized by Edney [43] is presented in reference [21]. For the convergence 
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study only one type of interaction, namely what is called the Type III interaction, is considered. 
Also the study is restricted to only one type of time stepping sequencing and only one value of the 
entropy correction parameter. The computational meshes (not shown) consist of 181 points in the 
circumferential direction and 91 points in the normal to the body direction and are highly clustered 
in the wall region to resolve the viscous layer. 

At this point, it is important to point out that the time step sequence used for the viscous steady 
flows is very different from the inviscid study. Most of the inviscid computations use the same initial 
time step input together with a local time-stepping procedure throughout the entire iteration process. 
The time step sequence chosen for the viscous steady calculations is based on experience with a wide 
range of hypersonic flow simulations and consists of doubling the time step every 100-400 time steps 
until the specified time step is reached. The initial time step is At = 0.001 which corresponds to a 
maximum Courant (CFL) number of 10 to 20 for the current problem and grid size. Larger values of 
the initial time step usually prevent convergence. The four specified time steps considered range from 
0.001 to 0.008 with the corresponding CFL numbers ranging from 20 to 200. Much larger maximum 
CFL (or specified time step) numbers are possible but do not improve the convergence rates. The 
value of the entropy correction parameter was fixed at S = 0.15, again based on experience with a 
wide range of hypersonic flow field simulations. 

The results of the blunt-body steady viscous flow obtained with the temporally second-order 
accurate algorithm (12) (hereafter referred to as the full matrix form ) are shown in figure 8. Here 
algorithm (12) for the viscous computations means the appropriate three-point central differences of 
the viscous terms are added to the explicit and implicit operators of (12). Depicted are the Mach 
contours ranging from 0 to 8 in increments of 0.1 and the entropy contours normalized with the 
freestream value and ranging from 0 to 6.4 in increments of 0.1. The final view in figure 8 is the L% 
norm residual history. The residual drops to machine accuracy (10~ 14 ) in less than 3200 steps. The 
corresponding results using the same algorithm (12a)-(12f) together with (13) (hereafter referred as 
the diagonal form) are illustrated in figure 9. No noticeable difference in the numerical resjfits is 
observed in the Mach number or entropy contours. However, the residual curves are very different. 
The residual for the diagonal scheme reached a plateau of 5 x 10~ 6 at 1500 steps and stayed at that 

level. 

A more difficult flow field computation is depicted in figures 10 and 11. The results using the 
same second-order time accurate full matrix algorithm are shown in figure 10. The convergence rate 
is slower than for the blunt body non-interfering case but is still satisfactory. The residual dropped 
seven orders of magnitude in 3000 steps. In both of the blunt-body flow with or without impinging 
shocks, steady-state can he reached within 1000-1500 iterations. The extra iterations are needed 
only to bring the residual to a lower level but no change in the coutour plots or surface pressures 
at least to within 3-4 digits of accuracy is observed. However, the results shown in figure 11 using 
the diagonal scheme are not satisfactory. The residuals dropped less than two orders of magnitude 
in 3000 steps. The noise appearing on the Mach number and entropy contours in the upper portion 
of the bow shock using the diagonal form of the scheme indicates that the algorithm has a problem 
reaching the converged steady-state solution. 

All of the results obtained for figures 8-11 have the viscous terms included in the implicit operator. 
If the viscous terms are not included in the impicit operator, then the full matrix scheme becomes 
unstable for At > 0.004, whereas the diagonal scheme exhibits no change in convergence rate. 

In summary, from the point of view based on the £ 2 -n° rni of the residuals, the best convergence 
rates were achieved by the full matrix form with the viscous terms included since it allowed the 
residual to drop to machine accuracy ( 10~ 14 ). The diagonal form (13) did not fare too well. Although 
there is a substantial savings in operation count per iteration, the Z 2 -norm of the residual never 
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dropped below 10 -6 for all the time steps considered. Moreover, the inclusion of the implicit viscous 
terms had little effect on the diagonal form of the scheme but is important for the full matrix form of 
the scheme. One way of taking advantage of the low operations count of the diagonal form (without 
the implicit viscous terms) is to use the scheme as an efficient way of obtaining a rough solution 
(from freestream) for the initial input to the full matrix algorithm. The temporally second-order 
time- differencing scheme had a marginal but beneficial effect on the convergence rates when compared 
with the temporally first-order scheme. 

Viscous Steady and Unsteady Mach 15 Computations With Impinging Shock Figure 12 illustrates 
the shock resolution of unsteady and steady thin-layer Navier- Stokes computations by the second- 
order time-accurate, full matrix algorithm (12). This steady test case is similar to the previous 
impinging shock study except the freestream Mach number is 15, the impingement shock angle 
is 22.75°, the freestream temperature is Too = 255. 6°K, and the Reynolds number based on the 
diameter is 186,000. Shown are the Mach contours from 0 to 15 in increments of 0.1. For the 
unsteady computation, the impingement shock at an angle of 22.75° relative to the freestream moves 
downward across the bow-shock of the blunt body. The impingement shock velocity is 10% of the 
freestream velocity (Moo = 15). Although the impingement shock locations for the unsteady and 
steady computations are similar, the shock patterns are very different. A 241 x 141 non-adaptive 
grid is used for both computations. A time step of 0.002 (equivalent to a maximum CFL of 48) is 
used for steady-state computations whereas a time step of 0.0005 (equivalent to a maximum CFL 
of 10-12 at the vicinity of the boundary layer and a CFL of 1 at the rest of the flow field) is used 
for the unsteady computations. The steady-state solution can be reached in 1200 steps with a three 
order of magnitude drop in the T 2 -norm residual. More detailed study of the surface pressure and 
heat transfer rate of these types of shock -on- shock steady and unsteady numerical simulations were 
reported in references [21,22]. 


VI. Concluding Remarks 

A two-parameter family of implicit time-accurate shock-capturing algorithms for hypersonic vis- 
cous flows has been developed. The proposed algorithms are formulated in finite volume and pseudo 
finite volume form and have been shown to be quite efficient and accurate for steady-state as well 
as unsteady viscous and inviscid hypersonic complex shock interactions. Some numerical aspects of 
these TVD-type algorithms that affect the convergence rate for hypersonic Mach numbers or real 
gas flows but have negligible effect on low Mach numbers or perfect gas flows are identified. Im- 
provements have been made to the algorithms to speed up the convergence rate in the hypersonic 
flow regime. Even with the improvements, though, the convergence is in general slower for real gases 
than for a perfect gas. The nonsmoothness in the curve fits of Srinivasan et al. may be a major 
contributing factor in slowing down the convergence rate. Also, the convergence rate is, in general, 
slower for viscous flows than for inviscid steady flows. 
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Fig. 1 The 31 x 33 grid. 
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Fig. 2 


The Mach contours (a), density contours (b), pressure contours (c) and k (d) computed by 
the implicit scheme (12) (6 = 1, u> = 0) for an equilibrium real gas with Moo = 15. 
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Fig. 3 


The Mach contours (a), density contours (b), pressure contours (c) and k (d) computed by 
the implicit scheme (12) (9 = 1, u> - 0) for an equilibrium real gas with M 0 0 = 25. 
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Fig. 5 Comparison of the I 2 -norm residual of a perfect gas and real gas computed by the scheme 
(12) (0 — 1, u? = 0) with Mqo = 25. Note that the scale of the ordinate for the perfect gas 
and the real gas is not the same. 
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Two-dimensional inviscid steady blunt-body flow with an impinging shock computed by the 
implicit scheme (12) (0 = 1, u> = 0) for a perfect gas with M = 4.6. 



29 





30 


Mach contours, entropy contours, and residual history for the steady viscous blunt body 
flow computed by algorithm (12,13) (6 = 1, u> = 1/2, diagonal form) with M = 8.03, 
Rco — 387,750, 7 - 1.4, and laminar flow. 





Fig. 11 Mach contours, entropy contours, and residual history for the steady viscous Type III 
shock interfernce flow computed by algorithm (12,13) (0 = 1, = 1/2, diagonal form) with 

Moo = 8.03, Rep = 387, 750, 7 = 1.4, and laminar flow. 
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